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ABSTRACT 



A bubble transport equation was developed and. a formal solution 
obtained for tbe bubble distribution- function.. The' distribution 
function obtained depends upon, the modet used, for; bubble' drag; .gas 
diffusion, and ocean circulation.. As a: specific application”, . solutions 
were obtained for a one-dimensional,. quiescent,, steady state ocean 
with a bottom source at 1 ST meters .. The: characteristic equations 
were integrated numerically to test” the: accuracy 'off approximate' 
analytical solutions.. The calculated bubble densities were compared 
with experimental data available for this problem.. When gas diffusion 
is neglected, the theoretical and: experimental distributions are in- 
general agreement- 
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I. niTRGDUCTrON 



Studies have been made by many* researchers regarding- the mechanisms 
of gas bubble dynamics in a fluid, from both the theoretical! [Ref . 1] 
and the experimental [Refs* 2, 3^ and. 4-1 aspects:.. The* mechanisms of 
gas bubble formation, transport y and extinction: are important: as. gas 
bubbles can effect the propagation of sound in the ocean, by the under- 
water scattering of sound beams.. Theoretical studies - have explained ' 
how a bubble varies in sine as. it rises in. a quiescent' fluid!. The 
phenomena involved includes shrinkage from gas diffusion, versus 
expansion from decreasing hydrostatic pressure and. shrinkage from 
surface tension for very small bubbles.. However,, there does not appear* 
to be at present any theory explaining how the distribution. of bubbles 
generated by a bottom or distributed, source in. the ocean, can be expected 
to vary with depth or initial bubble radius.. 

This paper outlines tbe theory for determining theoretically the 
distributions of bubbles generated, by a bottom source in the ocean 
[Ref. 5]. First, a linear bubble transport equation is developed and 
then it is solved by the method of characteristics.. A model for a one-- 
dimensional steady state ocean is specified for application* of this 
theory. Depending on the assumptions made regarding- the dynamics 
involved in the bubble motion,, such as the importance of gas diffusion 
or surface tension, various semi-analytical solutions, can.be obtained 
for the distribution of bubbles from a. bottom source.. These: distri- 
butions are then compared to available experimental data - .. In particular, . 
the distributions from a bottom source located at 197 meters were calcu- 
lated using various approximations and the results were compared with 
the experimental data in Ref 4.. 
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II. THE BUBBLE TRANSPORT EQUATION 



The transport equation is basically a balance equation [Refs. 6 
and 7]. It was desired to develop a transport: equation for the: 
description of the distribution of bubbles as a ftmctibn off position, 
velocity, size (bubble radius) and time in a circulating- field that' 
may contain sources or sinks* A distribution function that describes 
the bubble population can be defined and If the bubbles, satisfy the 
following general assumptions * this distribution: wil l , satisfy - a' 
Boltzmann-type transport equations 

(1) The bubbles are completely independent off one another.;. 

(2) If bubble interaction occurs r no more than, two hubbies will 
interact at any one time. 

Since bubbles in the ocean exact no force on one another and- the 
history of one bubble does not affect the history of another bubble, 
the two above assumptions hold and a transport, equation was found.. 

If i[>(r ,v,£ ,t )d 3 rd 3 vd£ is the number of bubbles at time t with 
radius dz about £, velocity d 3 v about v*,. and position d 3 r about r, 
then at time t f e t + dt, r r r y t f ) d 3 r r d 3 v f d£ T is the' number of 

bubbles at time t + dt with radius d£ f about £ r v velocity d 3 v f about v T , 
and position d 3 r f about r r , where r r y v T and z' are:- 

r T = r + vdt 

v r = v + adt (1) 

£* = £ t y)dt 

The above Taylor series expansions about time t neglects alll terms, of 
order two and higher. The changes in position, velocity, and radius 
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are the result of nonzero velocity v a dr/dtr, acceleration a =-dv/dt 5 . 
and radius change rate v = dj^/dt,. respectively.. 

Since the transport equation is a balance equation the number, of 
bubbles in a volume element d 7 x v =: d 3 r 'd^v'dX r at~ time: t: +- dt“, .must 
equal the number of bubbles that were in the volume: d^T' == d 3 rd 3 vd£ 
at time t plus or minus the bubbles that were generated or removed by 
any internal sources or sinks.. Thus an equation, for the: conservation* 
of bubbles is written as: 

\f;(r f ,v f ,£ f 5 t+dt)d 7 x r - ^(r r v r £ ,.t)d 7 x: = S(E*,v,jL,ir)d7r" (2) 

where S(r 9 v,£ ,t )d 7 x is the rate at which bubbles were introduced into 

m 

the volume element d'x about (r,.v 5 l) by distributed sources.. 

Expanding tJ; 1 = r r r v r r V ^t+dt) in. a Taylor series about.f ==i|j(r ) v } Jld) 
and keeping only terms through the first order yields:* 

ip* = $ + (3ij)/3t)dt -h (v + V ^)dt t (a. •* \MjOcLtr: +- v(dTj//9t)dt* (3) 

where 

*> 

V (•) 5 3CO/3X, a(-)/3y;: 3(~)/3z 
*> 

V ( ’ ) = 3C*)/3v , 3(-)/3v 3(-)/3v .. 

v x y z 

The volume element d 7 x r was transformed into terms of d 7 x by using 

a Jacobian. The results of this transformation yielded. 

7 = 3(P d 7 r (4 ) 

3( r r v t l) 

with 

> v ' »&' ) = 1 + ( v - a) dh +- ( 3v/3Jl)dt:.. (5) 

, -> -> •> ^ 

3 (r,v,£) 
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Substituting (3), (4) , and ( 5) into ( 2) dividing by d 7 xdt. 
and taking the limit dt 0 yields the bubble transport equation 
[Ref. 5] : 



3^ 

dt 



V ip + a 
r 



V ib 
v 



dip _ 
dt 



- <V J+ lr ) 



( 6 ) 



It should be noted that this transport: equation is: different: from 
the usual transport equation encountered because the term (V^-a _ t 3v/dt) 
is generally nonzero as the bubble radius change rate,- v(r, v,£ , t) , . and" 
acceleration, a(r,v,£,t), are functions odf bubble radius and' bubble 
velocity respectively. 

The distributed sources or sinks are contained, in the term 
S = S(r,v,£,t). This term can also include the rate at which bubbles 
are introduced into d t and removed from d t by scattering' off: small 
scale turbulent eddies or by bubble- bubble interactions.. Ib was assumed^ 
the S was independent of ip or proportional to ip.. This simpli cation 
yielded a first order, partial differential bubble transport equation . . 
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III. FORMAL SOLUTION TO THE. TRANSPORT EQUATION 



The bubble transport: Eq.. (6) irr the. form of- a. firnlr. order, partial 
differential equation was solved using the method of" characteristics 
[Ref. 5]. This was equivalent to replacing the transport: equation 
with a set of simultaneous y first order,, ordinary differential equations. 

The characteristic equations were determined- to: he the bubble 
dynamic equations: 



dr/dt = v 


r(t ) - v - (x. ,.y- ,z ) 
a a a o o: 


(7a) 


dv/dt = a 


v(t ) = v - (v ,v ,.v ) 

o o xo yo zo. 


(7b) 


di/dt = v 


z(t ) - Z 
a o 


(7c) 


d^/dt = S — 


0 

*3* 

ii 

/***> 

o 
,-P 
- — ✓ 


(7d) 



The model chosen determined the acceleration a(r,y,£,t) and rate of 
bubble radius change- v(r*v,.£ ,.t) thus establishing the. family of char- 
acteristic curves . This family was parameterized by the set off initial 

conditions (r ,v ) r and \b (r V v ) was the. initial distribution, 
o o o o o o o 

The general form of the distribution function was obtained by 
Integrating Eq. (7d) and making use again of a Jacobian.. The formal 
solution for the bubble distribution is- CRef .. 5]; 



X 



r,v,£,t)d 7 T = ^(r oT v o ,.£ o ,t) d 7 T Q t 



S(r fr 9 v Tr ,£ !l ,t J, )d 7 T ,T dt lt 



The integral here Is a line Integral along a characteristic, curve as 
are all the integrals In this thesis.. Equation (8) is: equivalent to 
saying that the number of bubbles in d7 T about - (r,v,£) atr. time: t - is. 
the number of bubbles that were originally in d 7 T Q plus or minus those 
added or subtracted by sources or sinks while proceeding along the 



characteristic curve . 
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IV. THE MODEL FOR THE SOLUTION OF THE BUBBLE TRANSPORT EQUATION 



Before solutions to Eq. (6) could be obtained, explicitly , . a - -. model 
of the upper ocean had to be specified.. This model, must provide-*- explicit 

expressions for the acceleration. a(r y v,l,.t) >• the rate of change of. 

— >■ — >■ — > — ->-■* 

bubble radius v(r,v,£,t), the bubble distribution S(r,v,£,t) and. the 

circulation field V(r,t). A Cartesian coordinate system was* chosen 

with the origin located at the surface of the ocean, and. the s-axis 

vertical and positive upward* Thus all depths were taken to* be. negative. 



A. THE BUBBLE RADIUS CHANGE RATE 

The rate of change of bubble radius for bubbles circulating in the 
ocean is caused by either gas diffusion a change irr the hydrostatic 
pressure, surface tension, or a combination of these, phenomena.. This 
interaction of effects was analyzed by assuming the bubbles to. be 
spherical containing an ideal gas* From the ideal gas* law 

(9) 



where R is the gas constant, T is the temperature,. P is the pressure, 
and V is the bubble volume. By taking the time derivative of'Eq. (9) 
and assuming isothermal changes, an expression for v 9 . the change in 
bubble radius was obtained: 

_ d£ = _ A AL 

v = dt 3 [n dt P dt 

The term — 4^ contains the gas diffusion effects while the effects of 
n dt 

the change in hydrostatic pressure and surface tension are.: contained 

. _ ^ 1 dP 

m the term ■=- . 

P at 



( 10 ) 
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The hydrostatic pressure P Is the sum of the' atmospheric - pressure , . 
P^, and the fluid pressure,, -p.^gz.. The surface tension pressure is 
2y/£ , where is the density of sea water,. p Q =1.03 x ID 3 ' Kg/meter 3 , 
and y Is the surface tension, of sea water,, approximately • equals to: 0.74- 
newt ons / me t e r . Therefore 



P - P - p gn + Zy/£ 
o o 



(ID 



and 



dP _ 2r d!L 

dt P o gV z i 2 - dt 



( 12 ) 



By using a simple model for gas diffusion [Ref. . 8] the following 
expression was obtained for dn/dt 

~ = - 64**3 ( P _ Fq ) (13) 

where 5 Is a constant. A value for 5 had. been determined from experi- 
mental data to be about 10“^ meters/second.. Dr was assumed for this 
model that the dissolved gas concentration gradient was confined to 
a thin shell at the air-water boundary because of the bubble f s' speed. 

An expression for v was obtained by substituting Eqs. (11) , . ( 12) , 
and (13) into Eq. (10) [Ref. 5j.. 



v = = l> l z +- 36RT'(z* - r)l/C'3(n-z)* + (14) 

at z 

where 

T = 2y/p Q g - 1.47 x 10 5 meters 2 

and 

D a P /□ g - 10 meters., 
o M o & 

If £ Is greater than 30 microns (1 micron (p) =10"^ meters), the 
surface tension terms can be ignored,. I.e. , f - 0. 



10 



B. THE BUBBLE EQUATION OF MOTION AND ACCELERATION 

A bubble that Is in motion with respect to a stationary reference 

system will experience several forces acting on it.. First, there is ' 

the force of the water acting on the bubble to pull it along. . Iff the 

bubble were completely entrained by the fluid so that the bubble's 

velocity (v) were equal to the circulation field velocity (V) the 

entrainment force would be p crdV/dt where cr it the volume of: the*-. 

o 

bubble . Complete entrainment can. occur' only* if? p: is: equal. to:p Q . . 

Since this is rarely the case* the bubble will maintain", a. relative 

— y ->• — ^ ^ ... 

velocity u = v - V with respect to the water.. This relative velocity 

— v m # 

produces a drag force, F , acting on the bubble to slow it. down*. 

As the bubble moves through the water v the bubble -pulls water 
along with It. The amount of water pulled along with it is equal, to 
some fraction, 6, of the bubble's volume.. From experimental data 3 is 
usually taken as 0.5. This phenomena is' called the mass adhering force 
and Is equal to -(Scrp^du/dt) .. The bubble will also experience a net 
bouyant force which is equal to crg(p Q — p) z v where g is the gravita- 
tional acceleration. 

By using Newton r s second law these forces were combined to: yield 

—V 

an expression for the acceleration a [Ref.. 5]: 



->■ 

F 



a = 



110 



dV 



(I-P/P ) 
o 






(S+p/pJ * (B+p/p ) dt * (B+p/p ) gZ 



(15) 



where 



dV SV vr ^ 

— = -r— 1 ( v-V )V •• 
dt &t r 



(16) 



In order to obtain an explicit expression for. the acceleration’ a, . 
it was necessary to derive an expression for the drag force F^. • Various 
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models for the drag force were considered, with the models being 
classified by using the Reynold r s number R^ = Zlup^/n where n is' 
the viscosity of water Cn 1Q~^ Newton— secorrd/meter^) . . 

For bubbles with a radius less than or equal', to 10D microns 
(£ _< 100 ), the Reynold r s number is less than one (R^ < 1) and:the 
drag is viscous according to Stokes Law:- 



\ a ^Stokes “ ( 

where K = 6^ in the presence of surface active material and in. the 
absence of these surface active materials.. 

The drag was considered to be both viscous - and turbulent - if : 

100p <_ 1500 u or 1 < < 7QQ- Thus the drag: force was 



-> 

F = F + F 

D Stokes Turb 



(18) 



Here ^ Sto j <:es was given as in Eq.. (17) with 6 tt~ < K < 12t~ and which 

arise from the energy dissipated in the turbulence: of the wake,, had the 
form: 



F Turb = -°- 5 Vo 5 ! uu (19) 

where is a drag coefficient and s^ is the bubble surface, area 

covered by the turbulent wake (s^ =' rri 2 -). If there were no surface 

. 

active materials present F_ would equal F„ , with K = 12 it. . 

* D H Stokes’ 

Since the oceans contain surface active materials the drag "force 
contains both a viscous drag and turbulent drag term. By noting that 
p/p^ << 1 and that a = 4tt£. 3 /3 for bubbles with i < 1500: y,. the: acceler- 
ation can be written as CRef. 53:: 



a = - (C-Vu - (£)uu + (^-)z + (1 + hdV/dt 
l 2 - l & B 



(20 
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with 



a = — — 1 8.75 x 10 6 meters /sec 2 

2 ep o 



and 



£ 3 K f I _ L 

5 = T “ft 7 = x. Id 1 

2 & 



If the bubbles were generated by a surface source such as rain, 
wave action, etc., the bubble radix are generally less than' 150 I y 
[Refs. 3 and 9] and the drag force is predominately- viscous . . Inr.this' 
case £ should be set equal to zero.. 

If the bubbles originate at a. bottom source, the radii are in** the 
range of 5Q0 y to 70Q y [Ref. 4] and the drag should be considered to be 
predominately turbulent with, a set equal to zero.. 



C. BUBBLE SOURCES 

The bubbles in the ocean might physically be considered to have 
originated from rhree types of sources.. The first type included 
bubbles that originated near the surface and. were caused by wave 
action, rain or snow. The distribution of bubbles originating"- from 
these phenomena was analyzed in Ref. 10.. 

The second type of bubbles are these nucleated by distributed 
particles of material as well as bubbles that originated, from biological 
activity such as photosythesis or bacteria and are distributed through- 
out the medium. 

The third type of bubbles includes those bubbles generated on the 
floor of the ocean. These bubbles are often sediment— initiated bubbles 
and measurements indicate these bubbles to he between 450 y — 800 y in 
radius with a peak around 550 y [Ref. 4].. It was this category of 
bubbles for which the distribution function was evaluated. It was 
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considered to be mathematically convenient to introduce the bubbles 
via boundary conditions rather than to include them in S(r ,v,£ 5 t) . . 
ip o ( r o V Q ^ & o ) represented the bubble distribution originating irr 
the bottom sediment. 

The model for the circulation field was chosen to be el quiescent 
fluid dictating that V = 0. For this model the bubbles^ could rise- 
only via buoyancy and the characteristic curves are vertical upward . 
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V. RESULTS 



For this thesis the bubble transport: theory was applied- to' a 
one-dimensional quiescent, steady state ocean' with a bottom source 
at 197 meters. The solution and approximate, solutions are given. . 

The characteristic equations were integrated, numerically to test" 
the accuracy of the approximate analytical solutions^. The calculated 
bubble distributions were compared with, experimental., data. . 

In this steady state model the velocities are vertical' upward, . 
and all functions depend spatxcally on. depth only.. The bubble transport 
equation (Eq. 6) can be written, as:* 



d6 d(t> 

v -r 1 - + a -y- 
dz dv 



dd 

v — --Ed 
d£ t^ 



t S_. • 



( 21 ) 



z 

Here <f.(z ,v,iL) is the distribution function,. S is the. distributed source 
function, v is the vertical velocity,, v is. the: rate of." change off bubble 
radius, and = 3a/ 3v t 3v/3£.. 

If z is chosen as the independent parameter,, the characteristic 
equations are [Ref. 5]: 



dv 

dz 



a_ 

v 



- \ v + f— 

£ gv 



d£ _ v _ £ 2 + 35RT( z£ - r)/v 
dz v 3( D-z) £ t 2 T 



and 



= _ 1 z a +Is.. 

dz v t Y v 



(22a) 

(22b) 



(22c) 



Since the characteristic curves are verti cal ly upward: along; the z- axis , . 

the characteristics or boundary conditions are v ~ v q ,. £ “ ^ Q >*a-nd 

= cf>(z ,v ,£ ) at z = z <0. 
o o o o — 
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Equation (22c) was integrated to yield 



z 



<Kz>v>£) 




‘o 



with the integration taking place along the characteristic curves 



Equation (24) is interpreted as a conservation equation which . states 
that the flux of bubbles, , at the point (z:,.v 5 l) is equal, to. the 
flux of bubbles, v dj , at the point (z V v v £ ) plus (or. minus.) any 
bubbles added (or substracted) by distributed sources along the 
characteristic curve. 

Prior to the calculation of the distribution function in either 
Eq. (23) or Eq. (24) it was necessary* to establish a mapping between, 
points (z,v,£) and (z r ,v T ,£ r ) on. the characteristic curve.. This was 
done by integrating simultaneously the characteristic equations (Eq. 
22) to obtain v(z,£; z r ,v r ,£ T ) and £( z ,£ ;z * ,.v r ,£ f ) for various 
z f 3 v f , and V and substituting these relationships into Eq. (23). 
Since ct/£ 2 >> 1, it was very time consuming to integrate. Eq.. (22a.) 
numerically. Referring to Fig- 1^ it is seen that an. appropriate 
approximation is to equate the velocity and the terminal velocity. . 
This permits semi-analytical solutions.. 



{z,v(z;z ,v ,£ ),£(z;z ,v ,£ )*, z < 0}.. By- using- the Jacobian*. Eq. (23) 
o o o o o o — 

was rewritten as [Ref. 51: 



vd(z ,v,£ )dvd£ = v dj ( z r v ,£ )dv* d£ 




(24) 
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This approximation can be verified by rewriting Eq. (22a) as 
[Ref. 5 3 



vrhere 



and 



- — — ( — ) ( v — v ) ( v -}- v 1 ) 

dz T T. 



V T 



-1 + Cl + Hv t1 /v T2 ) 2 ‘3 



V = 



2 ( v /v ^ ) 
T1 T2 



1 4- [1 + 4(v T ^/v T2 ) 2 -J 



Y 2 - 



1/2 



2(V TL /v T2 2) 



( 25 ) 



(26) 



(27) 



T1 



It 

8 a 



(28a) 



is the terminal velocity in. the case where viscous' drag dominates 
(£ < 100 u ). 




(28b) 



is the terminal velocity in the case where turbulentr drag dominates. 
(jL > 400 ji). Thus 







for 

for 



Z £ LOO 
l > 400 



(29) 



Equation (25) was integrated using fourth rank optional. Runge Kutta 
[Ref. 11], These numerical results are presented graphically in Fig. 1. 
It is observed that when a bubble is released' from a depth of- 20.00 
meters with zero velocity y it will reach terminal velocity after having 
traveled a distance of only several bubble radii'.. Thus the approximation 
that v = v^ is considered to be a valid approximation for this model. 

This approximation is convenient as Eq. (22a) is very difficult to 
integrate numerically. 
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By using the appro ximat ion v = v^ y the bubble distribution function 



can be written approximately as r. 



where 



4)(z,v,£) = <5(v - v )$( z v £) 



$(z,£) = / c(i(z y v y £)dv 



(30) 



(31) 



is the "depth-radius" distribution function.. By substituting- Eq. . ( 31) 
into Eq. (24) and integrating- over v 5 . an expression for $ was determined 
as [Ref. 5] : 

( ar 



z 






e / dZ ’ <32) 



where 



and 



s(z* ,2, ' ) = 



dl 



-m 

o \3£ / 



S( z r y v r y £ r )dv r 



d£ y etc.. 



(33a) 



(33b) 



In order to calculate $(z y £) y it was necessary to~ determine 3 £ T /g£ 
and £ f (z' ;z,£) either analytically or numerically. Note that 
dZ'/d& = l/(8£/3£ r ), and that the approximate differential equation 
that 8y3£’ satisfies is obtained by setting v = in Eq.. (22b) and 
taking the partial derivative with, respect to £ r [Ref.. 5] : 

r 



a a£ 
a z 3£ 



m [-( 

H™)H 



* i ■ ■ 



fd V 



z Ui —I 



/ [«»-> 

/ 



3(U-z) -t- ~ 



( 34a) 



where 



-r/ 2 



3v t _ 3v ti 

32. 2. 



1+4 



v V 2 
T1 

J 

T2 



~r •• 



(34b) 



Furthermore, with the approximation v = , then 2.= 2.^, and Eq. (22b) 
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becomes 



(35) 



l 2 + 3<5RT( zl - D/v 



dz 



3(D-z U +~ 2T 



Equations (34) and (35) can be Integrated, numerically using optimal 
fourth rank Runge Kutta [Ref. 1X3- These results can be substituted 



a good approximation. It Is expected that the approximation l -- will 
also be very good. 

To obtain approximate analytical results for comparison, to - the 
more exact numerical results outlined above,, the approximations of: 

Eq. (29) can be used. Neglecting surface tension,. TT =- 0 is also a. 
valid approximation for £ >_ 30 u .. Then Eq.. (31) can be integrated 
analytically [Ref. 53. For v - v and r - 0 



into Eq. (32) to obtain $ = for various sources.. Since v=-v^ is 




1/3 



(36a) 



L 



For v = v^ 2 and F = 0 





2/3 



(37a) 



The inverse transformation are 




(36b) 



and 




(37b) 



respectively. Then 



32.' _ ,D-z . I 1 

"51 H 



(36c) 
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3£ f 
3 £ 



(37c) 




respectively. If gas diffusion, can be neglected,, 6 = 0,then Eqs..(36) 
and ( 37 ) reduce to : 



_£ 

V 



3 £* 
3 £ 




1/3 



1/3 



(38a) 

(38b) 



Figures (2)-, (3), (4) T (5), and (6) demonstrate the results obtained 
by integrating Eq. (35) numerically and comparing £^ to the analytical 
approximations of Eqs. (36), (37) ^ and. (38) .. Figures (2) , . ( 3) , . and (4) 
were drawn for a depth of 20 meters and Figs.. (5) and (6) were drawn' 
for a depth of 197 meters* Figures (2) and (50 demonstrate* that £^ 
is a good approximation for £^ when £ <r 100 p,.. However,, this approxi-- 
mation breaks down when £ < 30 p-.. This is because the effects of 
surface tension, which are important in bubbles of~ £ <__30 p, .are 
neglected. Figures (3), (4), and (6) demonstrate that £^2 is a good 
approximation for £^ when £ >_ 400 p.. From Fig.. (6),. it is. observed 
that £ is not a good approximation for £.. This is because £ 

I J i o 

neglects the effects of gas diffusion. 

The survival threshold for bubbles Is different depending on 
whether the bubbles were released at a depth 20 meters or 197 meters. 

In the 20 meter case, bubbles of £j> LOO p will reach the. surface , . 
while in the 197 meter case, bubbles of £>^1000 meters will. surface. . 
This difference Is caused by gas diffusion being dominate over: the 
decrease in hydrostatic pressure- Since gas diffusion is proportional 
to (P-P )/P, it will be constant for z < 30 meters and will only 



decrease for z _> 30 meters. Thus bubbles released a"t 197 meters are 
exposed to a high gas diffusion rate for practically alll their travel, 
while bubbles released at 20 meters are exposed to a decreasing gas 
diffusion rate. 

Finally the bubble distribution, function was - calculated- for each 
of these approximate solutions to Eqs.. (22a) and (22b), (3M-) and (35),. 
(36), (37), and (38). It was convenient to introduce the initial- 
distribution function as the boundary condition. z;zr^ ,£ q ) instead- 
of as a distributed source. Since s(n r ,Jt r ) is zero,, the solution can 
be written as 



7 



* (z > z o’V = v r 



To 



*(z r i ) i/On/an'). 
o o 



(39a) 



For v = v , Eqs. (34) and (35) must be integrated numerically to yield 



$_(z\z ,£ ). For v =v (2 < IQQu) and E - 0 

I 3 O O II — 



Ti 



$ - $ 



for v ^Q^i) 



(z;z ,2. ) - (-tf) $(z ,X );■ 
1 o o ^D— z^y o o 



4 = y^vv =(& 



1/2 



<£ (z ,1 ); 
o a 



(39b) 



(39c) 



and if the effects of gas diffusion are not considered 

V I \ 1/3 

$ = $ „(z;z ,i ) = tt^$(z ,-Z ) ( u — ~ — ) 

3 o o V m o o Vb-z / 



(39d) 



Here V is computed using 

1 O 

£ = £ 



D-z 



1/3 



o D-z 



in Eq. (26). 

The above distributions were calculated using an initial distribution 



function, $(z ,£ ), obtained from Ref. 4 . The authors of Ref 1 4 
o o 
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determined the distribution of bubbles rising from "the highly organic, 
anaerobic sediment on the bottom 11 of Saanich Inlet,. British Columbia. 

The depth of the water was 197 meters.. The density of these rising 
bubbles was determined by using stationary echo-sounders operating 
at various frequencies in the kilohertz' range.. Figure 7a is a' histogram* 
that presents the presumed bubble distribution at ai depth of 
197 meters as presented in Ref* 4* Figures (7b) through (7g) demon- 
strate how this initial distribution changes' as the bubbles' rise in- 
height. It should be noted that height is measured from the ocean 
bottom upward to the surface. 

The different histograms at each height above the. ocean floor shows 
how the bubble distribution is changing depending on the; model chosen 
to be used for the calculation of £^.. In the model where gas' diffusion 
as well as hydrostatic pressure changes are considered^,. the bubble 
distribution $(z;z o ,.£ q ) shifts tpward smaller size bubbles as. height 
increases. It can also be observed that the width of the intervals 
is increasing with height. Finally all the bubbles disappear entirely 
at a height of about 100 meters. 

In the small bubble approximation $ - z;:z^ ,J^) ,. where the drag 

force is considered to be viscous and v = v pq 5 * the distribution 

appears to spread out by shifting toward bubbles of both smaller, and 

larger radii. Bubbles with initial radius less than 900 y continually 

decrease in radius, while bubbles with an initial radius greater than 

900 y will initially decrease in size but as the effects of gas 

diffusion reduce, these bubbles expand.. In the turbulent, drag model, . 

v = v__ , $ = $_(z:z , £ ) shifts toward smaller radius bubbles because 
i Z Z o o 

of the effect of gas diffusion and disappears at a height of about 125 
meters. 
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When gas diffusion is neglected, the bubble distribution 

$ * (z:z ) shifts toward bubbles of larger radii.. This was to 
3 o o 

be expected because as hydrostatic pressure decreased,, the bubbles 
would increase in size. The bubble dislzribution caicnlatredl for 
^^(z;z^,^ o ) compares quite well with the experimental data of 
McCartney and Bary in Ref. 4. It appears as though gas diffusion 
was somehow not important in the bubbles they measured.. These 
bubbles originated at the bottom, so perhaps they had some k-ihd of 
coating that impeded gas diffusion across their surfaces.. 
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VI. CONCLUSION 



Calculating bubble distributions in the ocean by solving- the 
bubble transport equation appeared to be quite feasible.. Although 
the solution of the transport equation can be quite iirrolvecL,. a careful 
choice of approximations often permits approximate analytical solutions . 
Otherwise, the equations can usually be numerically' integrated.. For 
the model chosen, analytical expressions were obtained for several 
approximations. All of these expressions yielded, results that: were 
in general agreement with experimental data or with intuitive reasoning.. 
It was seen that the drag force model chosen and the inclusion of gas 
diffusion greatly influenced the results obtained.. 

The assumption that the velocity of the bubble was equal to the 
terminal velocity should be quite valid as the data presented on Fig.. 

1 demonstrated. With this approximation, one can. easily calculate the 
bubble radius at any depth. The radius is ■ strongly dependent on the 
model chosen for the drag force, gas diffusion and surface tension.. 

The choice of the drag model determines the equation used to calculate 
the terminal velocity. 

Finally the distribution function for a bottom source was calculated. 
The data generated by this theoretical approach agrees well with the 
limited experimental data available. The distribution calculated for 
the no gas diffusion model yields data that is in good agreement with 
the experimental data in Ref. 4. 

From the limited amount of experimental data available,, it appears 
that gas diffusion is unimportant in bubbles originating in bottom 
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sediment and transported to the surface- It is felt that:, the type of 
gas in the bubbles would greatly influence the amount of diffusion that 
would take place [Ref. 12]. Since a fluid* super-saturated with the 
same gas that is within the bubble, is the only fluid, in which, gas; 
diffusion is not considered the dominating effect,, there mustr he_ some 
mechanism inhibiting gas diffusion from sedimentr bubbles .. Further 
investigation into this phenomena seems warranted. 
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HEIGHT OFF BOTTOM IN METERS 




CHANGE IN BUBBLE RADIUS AS A FUNCTION OF DISTANCE 
BUBBLE HAS RISEN FROM 2,=- 197 METERS A COMPARISON 
OF APPROXIMATE ANALYTICAL SOLUTION TO EXACT 
NUMERICAL SOLUTION. 



Figure 5 
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BUBBLE RADIUS IN MICRONS (aO 



5000- 



L = 

LT,= 



lt 2 = 




CHANGE IN BUBBLE RADIUS AS A FUNCTION OF DISTANCE 
BUBBLE HAS RISEN FROM 2 0 =-l 97 METERS A COMPARISON 
OF APPROXIMATE ANALYTICAL SOLUTION TO EXACT 
NUMERICAL SOLUTION. 



Figure 6 



NUMBER OF BUBBLES 



DISTRIBUTION OF BUBBLES ORIGINATING IN 



SEDIMENT AT DEPTH 197 METERS. 




INITIAL BUBBLE DISTRIBUTION i ( E,L) (REF. 4) 



Figure 7a 
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NUMBER OF BUBBLES NUMBER OF BUBBLES 





I 2 = 

h' 



BUBBLE DISTRIBUTION AT HEIGHT = 25 METERS . 



F igure 7b 
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BUBBLE RADIUS (>u) 




*r 

X 2 = 

£3= 



BUBBLE DISTRIBUTION AT HEIGHT = 50 METERS. 



Figure 7c 
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BUBBLE DISTRIBUTION AT HEIGHT = 75 METERS. 
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Figure 7d 
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BUBBLE DISTRIBUTION AT HEIGHT = 100 METERS 



Figure 7e 
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Figure 7f 




BUBBLE DISTRIBUTION AT HEIGHT = 197 METERS . 



Figure 7g 
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COMPUTER OUTPUT 



HEIGHT OFF THE BOTTOM 



XLO 


RADIUS 


PH 


400.000 


0.0 


0. 


465.000 


240.062 


0. 


475.000 


253.600 


0. 


485.000 


276. 132 


0. 


490.000 


284.593 


6. 


500.000 


301.010 


6. 


510. 000 


316. 856 


6 . 


515. OOC 


324.596 


31. 


530.000 


347.201 


31. 


545.000 


369. 032 


31. 


555.000 


383.230 


15. 


575.000 


410.924 


15. 


595.000 


437. 342 


15. 


600.000 


444.469 


3. 


635.000 


489.908 


3 . 


660.000 


521. 523 


3. 


670.000 


534.007 


0. 


705. 000 


577. 089 


0. 


730.000 


607.364 


0. 


740.000 


619.375 


0 . 


775.000 


661.022 


0. 


800.000 


690.444 


0. 


810.000 


702.146 


0. 



XLO 


RLT1 


PH 


400.000 


338.683 


0. 


465.000 


431.628 


0. 


475.000 


444.714 


0. 


485.000 


457.588 


0. 


490.000 


463. 954 


6. 


500.000 


476.556 


6. 


510.000 


439.000 


6. 


515.000 


495. 169 


29. 


530.000 


513.480 


29. 


545. 000 


531.533 


29. 


555.000 


543.442 


14. 


575 .000 


566.997 


14. 


595.000 


590.247 


14. 


60C.C00 


596.018 


3. 


635.000 


636.015 


3 . 


660. 000 


664.219 


3 . 


670.000 


675.430 


0. 


705 .000 


714.391 


0. 


730. 000 


741.996 


0. 


740.000 


752.992 


0. 


775.000 


791.303 


0. 


800.000 


318.520 


0. 


810. OCO 


829.376 


0. 



Z5.QGQ METERS 



RLT3 


FHI3L 


417.535 


0.0 


485.385 


0 ..929 


495.823 


0..929 


506.261 


0.930 


511.480 


6. .509 


521.919 


6.510 


532.357 


6.512 


537.576 


3L.633 


553.234 


31.643 


568.892 


3L.653 


579.330 


15.830 


600.207 


15.835 


621.083 


13 ..840 


626.303 


3.727 


662.837 


3:. <729 


688.933 


3. .730 


699.371 


0. 933 


735.906 


0.933 


762. 002 


0.933 


772.440 


0.933- 


808.974 


0.934 


835. 070 


0.934 


845.509 


0.0 



RLT2 


phi:z 


206.200 


0.0 


297. 174 


0 ..93 8 


310. 439 


Q.93 8 


323.562 


0.-938 


330.075 


6.564 


343. 009 


6 .564 


355.830 


6.564 


362.201 


31.881 


331.168 


31.881 


399.935 


31.881 


412.346 


15.940 


436. 950 


15.940 


461.293 


15.940 


467.342 


3.751 


509.318 


3.751 


538.952 


3.751 


550. 735 


0.938 


591.636 


0.938 


62C.693 


0.938 


632.244 


Q. 938 


672.467 


0.938 


701.018 


Q.938 


712.401 


0.0 



I 

0 

9 24 

925 

926 

482 

486 

490 

530 

552 

572 

792 

802 

812 

721 

724 

726 

932 

932 

933 

933 

933 

933 

0 

II 

0 

879 

879 

879 

155 

155 

155 

394 

894 

894 

947 

947 

947 

517 

517 

517 

879 

379 

879 

879 

879 

379 

0 
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HEIGHT OFF THE BOTTOM = 50.000 METERS 



XLO 

400.000 

465.000 

475.000 

485.000 

490.000 

500.000 

510.000 

515.000 

530. 000 

545.000 
555 .000 

575.000 

595.000 

600.000 

635.000 

660.000 

670.000 

705.000 

730.000 

740.000 

775.000 
800 .COO 

810.000 



XLO 

400.000 

465.000 

475.000 

485.000 

490.000 

500.000 

510. 000 

515.000 

530.000 

545.000 

555.000 

575.000 

595.000 

600.000 

635.000 

660.000 

670.000 

705.000 

730.000 

740.000 

775.000 

800.000 

810.000 



RADIUS 
0 . 0 
0.0 
0.0 
0 . 0 
0.0 
0 . 0 
0. 0 
0.0 
3 . 0 
0. 0 
0.0 

92.385 
188.945 
204.448 
290.264 
340.102 
358.626 
419. 267 
459.715 
475.396 
528. 524 
565.120 
579.499 



RLT 1 
222.533 
389. 746 

408.047 
425.532 
434.014 
450.529 
466. 521 
474.345 
497.222 
519. 336 
533.722 
561.775 
589. 029 
595.736 
641.693 
673. 641 
686.252 
729.760 
760.329 
772.455 
814. 509 
844.226 

856.048 



PHI 
0. 0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 

14.316 

14.403 
3.393 
3.410 
3.419 
0.855 
0. 857 
0.858 
0.359 
0.860 
0.361 
0.0 



PHI1 
0.0 
0.758 
0.758 
0.758 
5.309 
5.309 
5.3 09 
25.787 
25 .737 
25.787 
12.894 
12.894 
12.894 
3.034 
3.034 
3.034 
0.758 
0.758 
0.758 
0.758 
0.758 
0.758 
0.0 



RLT3 
438.615 
509. 890 
520.855 
531. 821 
537.303 
548.269 
559.234 
564. 717 
581.165 
597. 613 
608. 578 
630. 509 
652.440 
657.923 
696.301 
723. 715 
734.680 
773. 059 
800.473 
811.438 
849. 817 
877.230 
888.196 



RLT2 
0.0 
0. 0 

48. 238 
79.764 
93. 391 
118.240 
140.916 
151.662 
182.147 
210.647 
228. 815 
263.610 
296. 797 
304. 888 
359.687 
397. 238 
411.959 
462.363 
497.464 
511.326 
559.139 
592. 707 
606.015 



PHI 3 
0.0 
0.855 
0.855 
0.856 
5.991 
5.994 
5.997 
29. 136 
29.155 
29.174 
14. 593 
14.603 
14.613 
3.439 
3.443 
3.445 
0.861 
0.862 
0.863 
0.863 
0.863 
0. 864 
0.0 



PH I 2 
0.0 
0.0 
0.871 
0.871 
6.096 
6. 096 
6.096 
29.610 
29.610 
29.610 
14. 805 
14. 805 
14. 805 
3.484 
3.484 
3.484 
0. 871 
0.871 
0. 871 
0. 871 
0.871 
0. 871 
0.0 
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HEIGHT 


OFF THE 


BOTTOM = 


75.000 


METERS 


XL 0 


RADIUS 


PHI 


RLT3 


PHI 3 


400.000 


0.0 


0.0 


464. 721 


0.0 


465.000 


0.0 


0.0 


540.238 


0.775 


475.000 


0 . 0 


0.0 


551.856 


0.776 


485.000 


0.0 


0.0 


563. 474 


0.777 


490 i 000 


0.0 


0.0 


569.283 


5.440 


500. 000 


0 . 0 


0.0 


580. 901 


5 . 444 


510.000 


0.0 


0.0 


592. 519 


5.449 


515.000 


0.0 


0.0 


598.328 


26.474 


530. 000 


0 . 0 


0 . 0 


615. 755 


26.502 


545.000 


0.0 


0.0 


633.182 


26. 528 


555.000 


0 . 0 


0.0 


644.800 


13.272 


575. OCO 


0 . 0 


0.0 


668. 036 


13.288 


595.000 


0.0 


0.0 


691.272 


13.302 


600. 000 


0 . 0 


0.0 


697.081 


3.131 


635.000 


0.0 


0.0 


737. 744 


3. 136 


660.000 


0.0 


0.0 


766. 789 


3. 139 


670. 000 


0 . 0 


0.0 


778.407 


0.785 


705.000 


148.353 


0.766 


819. 070 


0.786 


730.000 


244.821 


0.772 


848. 115 


0.787 


740. COO 


273.440 


0.774 


859. 734 


0.787 


775.000 


358.003 


0 . 778 


900. 397 


0.788 


800. 000 


410.244 


0.780 


929.442 


0.788 


810.000 


429. 964 


0.0 


941.060 


0.0 



XLO 


RLT1 


PHI1 


RLT2 


P H I 2 


400. COO 


o. o 


0.0 


0 . 0 


0.0 


465.000 


331.106 


0.638 


0.0 


0.0 


475.000 


360.093 


0.638 


0.0 


0.0 


485. OCO 


386. 043 


0.638 


0 . 0 


0.0 


490.000 


398.163 


4.464 


0 . 0 


0.0 


500.000 


421.059 


4.464 


0.0 


0.0 


510.000 


442.507 


4.464 


0 . 0 


0.0 


515.000 


452.784 


21.681 


0.0 


0.0 


530. 000 


482.162 


21 .681 


0 . 0 


0.0 


545.000 


509.784 


21.681 


0 . 0 


0.0 


555.000 


527.421 


10.841 


0.0 


0.0 


575.000 


561.206 


10. 841 


0 . 0 


0.0 


595.000 


593.406 


10. 841 


12. 645 


13.575 


600.000 


601.252 


2.551 


42.242 


3.194 


635.000 


654.335 


2. 551 


153.845 


3.194 


660.000 


690.674 


2.551 


212. 749 


3. 194 


670.000 


704.917 


0.638 


234.349 


0.799 


705 .000 


753. 684 


0.638 


304.416 


0.799 


730.000 


787.663 


0.638 


350.758 


0.799 


740. 000 


801.088 


0.638 


368.659 


0.799 


775.000 


847.439 


0.638 


429. 073 


0.799 


800.000 


880.028 


0.638 


470. 501 


0. 799 


810. 000 


892.959 


0.0 


486.742 


0.0 
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HEIGHT OFF THE BOTTOM = 100.000 METERS 



XL 0 

400. COO 

465.000 
475. OOC 

485 .000 

490. 000 

500.000 

510.000 

515.000 

530.000 

545.000 

555.000 

575.000 

595.000 

600.000 

635.000 

660.000 

670.000 

705.000 
730. 000 

740.000 

775 .000 
800. 000 

810.000 



XLO 

400.000 

465.000 

475.000 

485.000 

490.000 

500. 000 

510.000 

515.000 

530.000 

545.000 

555.000 

575.000 

595.000 

600. 000 

635.000 
660. COO 

670.000 

705.000 
730. COO 

740. 000 

775.000 

800.000 

810.000 



RADIUS 
0.0 
0 . 0 
0.0 
0.0 
0.0 
0.0 
0.0 
0 . 0 
0.0 
0.0 
0. 0 
0.0 
0.0 
0. 0 
0.0 
0. 0 
0. 0 
0.0 
0.0 
0.0 
0.0 

136.578 

194.535 



RLT1 

0.0 

219.089 

285.770 
332.348 
351.978 
386.703 
417.224 
431.336 
470. 296 
505.502 
527.438 
568. 553 
606.878 
616.113 

677.771 
719.318 
735.487 
790.452 
828.450 
843.406 
894.835 
930.829 
945. 079 



PHI 
0.0 
0.0 
0.0 
0.0 
0.0 
0. 0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
C. 0 
0.0 
0.0 
0.0 
0.0 
0.680 
0. 0 



PH 1 1 
0.0 
0.517 
0.517 
0.517 
3.618 
3.618 
3.618 
17. 575 
17. 575 
17.575 
8.787 
8.787 
8.787 
2.068 
2.068 
2.068 
0.517 
0.517 
0.517 
0.517 
0.517 
0.517 

o.o- 



RLT3 
498.412 
579.404 
591. 865 
604.325 
610.555 
623. 016 
635.476 
641. 706 
660. 396 
679.087 
691.547 
716.468 
741.388 
747.619 
791.230 
822.380 
834. 841 
878. 452 
909.603 
922. 063 
965. 674 
996.825 
1009. 285 



RLT2 
0. 0 
0. 0 
0.0 
0 . 0 
0.0 
0 . 0 
0 . 0 
0.0 
0 . 0 
0.0 
0.0 
0 . 0 
0 . 0 
0.0 
0 . 0 
0 . 0 
0.0 

54. 319 
147.833 
177.412 
267. 406 
324. 165 
345.729 



PHI 3 
0. 0 
0.690 
0. 691 
0.692 
4.846 
4. 851 
4. 856 
23.599 
23.634 
23.666 
11.843 
11.862 
11.880 
2.796 
2. 803 
2.807 
0.702 
0. 703 
0.7 04 
0.704 
0. 705 
0.706 
0. 0 



PHI 2 
0. 0 
0. 0 
0.0 
0. 0 
0.0 
0.0 
0.0 
0. 0 
0.0 
0. 0 
0.0 
0.0 
0.0 
0.0 
0. 0 
0. 0 
0. 0 
0. 719 
0.719 
0.719 
0.719 
0.719 
0.0 



41 



HEIGHT 


OFF THE 


BOTTOM = 


125.000 


METERS 


XL 0 


RADIUS 


PHI 


RLT 3 


PH I 3 


400.000 


0.0 


0.0 


544.643 


0.0 


465.000 


0 . 0 


0.0 


633. 148 


0.596 


475.000 


0.0 


0.0 


646 . 764 


0.597 


485.000 


0.0 


0.0 


660.380 


0.598 


490.000 


0 . 0 


0.0 


667.188 


4.191 


500.000 


0.0 


0.0 


680. 804 


4 . 197 


510. 000 


0.0 


0.0 


694.420 


4.203 


515.000 


0 . 0 


0.0 


701.228 


20.428 


530.000 


0.0 


0.0 


721.652 


20.467 


545.000 


0.0 


0.0 


742.076 


20.504 


555.000 


0 . 0 


0.0 


755.693 


10.264 


575.000 


0.0 


0.0 


782.925 


10.286 


595.000 


0 . 0 


0.0 


810. 157 


10.306 


600.000 


0.0 


0.0 


816.965 


2.426 


635 .000 


0.0 


0.0 


864. 621 


2.433 


660. 000 


0.0 


0.0 


898.661 


2.438 


670.000 


0.0 


0.0 


912.278 


0.610 


705.000 


0.0 


0.0 


959.934 


0.611 


730.000 


0 . 0 


0 . 0 


993. 974 


0.612 


740.000 


0.0 


0.0 


1007 . 590 


0.613 


775.000 


0.0 


0.0 


1055.246 


0.614 


800.000 


0 . 0 


0.0 


1089.287 


0.614 


810.000 


0.0 


0.0 


1102.903 


0 . 0 



XL 0 


RLT 1 


PHU 


RLT 2 


PH I 2 


400.000 


0.0 


0.0 


0.0 


0.0 


465.000 


0.0 


0.0 


0.0 


0.0 


475 .000 


0 . 0 


0.0 


0 . 0 


0.0 


485.000 


237.512 


0.396 


0.0 


0 . 0 


490.000 


281.882 


2.773 


0.0 


0.0 


500.000 


344. 694 


2.773 


0 . 0 


0.0 


510.000 


392.070 


2.773 


0.0 


0 . 0 


515.000 


412.549 


13 .469 


0.0 


0.0 


530 . GOO 


466.053 


13.469 


0 . 0 


0 . 0 


545.000 


511.787 


13.469 


0.0 


0 . 0 


555.000 


539.423 


6.734 


0.0 


0.0 


575.000 


589.940 


6.734 


0 . 0 


0 . 0 


595.000 


635.914 


6.734 


0.0 


0.0 


600.000 


646. 868 


1.585 


0 . 0 


0.0 


635.000 


719.049 


1.585 


0 . 0 


0 . 0 


660.000 


766.966 


1.585 


0.0 


0.0 


670.000 


765. 493 


0.396 


0 . 0 


0.0 


705.000 


848.073 


0.396 


0 . 0 


0 . 0 


730 . OCO 


891.037 


0.396 


0.0 


0.0 


740.000 


907. 893 


0.396 


0 . 0 


0.0 


775.000 


965.657 


0.396 


0.0 


0 . 0 


800.000 


1005.929 


0.396 


1 14.088 


0.629 


810.000 


1021.843 


0.0 


153. 122 


0 . 0 



42 



HEIGHT OFF THE BOTTOM 



XLO 


RADIUS 


PH 


400.000 


0.0 


0 . 


465.000 


0.0 


0 . 


475 .000 


0.0 


0 . 


485.000 


0.0 


0 . 


490.000 


0.0 


0 . 


500.000 


0.0 


0 . 


510 . 000 


0.0 


0 . 


515.000 


0 . 0 


0 . 


530.000 


0.0 


0 . 


545.000 


0.0 


0 . 


555.000 


0.0 


0 . 


575.000 


0.0 


0 . 


595. 000 


0 . 0 


0 . 


600.000 


0.0 


0 . 


635.000 


0.0 


0 . 


660.000 


0 . 0 


0 . 


670.000 


0.0 


0 . 


705.000 


0.0 


0 . 


730.000 


0 . 0 


0 . 


740.000 


0.0 


0 . 


775. 000 


0 . o • 


0 . 


800.000 


0 . 0 


0 . 


810.000 


0.0 


0 . 



XLO 


RLT 1 


PH 


400.000 


0.0 


0 . 


465. 000 


0.0 


0 . 


475.000 


0 . 0 


0 . 


485.000 


0.0 


0 . 


490. 000 


0.0 


0 . 


500.000 


288.314 


1 . 


510.000 


372.655 


1 . 


515. 000 


404. 24 5 


9 . 


530.000 


480. 123 


9 . 


545.000 


540.491 


9 . 


555.000 


575.747 


4 . 


575.000 


638.591 


4 . 


595.000 


694.510 


4 . 


600.000 


707. 697 


1 . 


635.000 


793.635 


1 . 


660.000 


349.981 


1 . 


670.000 


871. 653 


0 . 


705.000 


944.486 


0 . 


730. 000 


994.222 


0 . 


740.000 


1013.686 


0 . 


775.000 


1080.214 


0 . 


800. 000 


1 126. 462 


0 . 


810.000 


1144.711 


0 . 



150.000 METERS 



RLT 3 


PH I 3 


614.835 


0.0 


714. 746 


0.490 


730 . 116 


0.491 


745.487 


0.492 


753. 173 


3.445 


768. 544 


3.451 


783.914 


3.458 


791. 600 


16.809 


814. 656 


16. 851 


837. 713 


16.890 


853. 083 


8.457 


883.825 


8.481 


914 . 567 


8.502 


922.252 


2.002 


976.050 


2.009 


1014.478 


2.014 


1029. 848 


0 . 504 


1083.646 


0.506 


1122. 074 


0.506 


1137.445 


0 . 507 


1 191.243 


0.503 


1229.670 


0.509 


1245.041 


0 . 0 



RLT 2 


PH I 2 


0.0 


0 . 0 


0.0 


0.0 


0 . 0 


0.0 


0.0 


0.0 


0.0 


0.0 


0 . 0 


0.0 


0 . 0 


0.0 


0.0 


0.0 


0 . 0 


0 . 0 


0.0 


0.0 


0 . 0 


0.0 


0 . 0 


0.0 


0.0 


0.0 


0 . 0 


0.0 


0.0 


0 . 0 


0.0 


0.0 


0 . 0 


0.0 


0.0 


0 . 0 


0.0 


0.0 


0 . 0 


0 . 0 


0.0 


0.0 


0.0 


0.0 


0 . 0 


0.0 



II 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 1 

0 

0 

0 

0 

0 

928 

928 

362 

362 

362 

681 

681 

681 

101 

101 

101 

275 

275 

275 

275 

275 

275 

0 
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HEIGHT 


OFF THE 


BOTTOM = 


E75..000 


METERS ' 


XLO 


RADIUS 


PHI 


RLT3 


PHI!3 


400. 000 


0.0 


0.0 


745 . 304 - 


O'.O 


465. GOO 


0.0 


0.0 


866.416 


0..360 


475.000 


0.0 


0.0 


885.048 


0 ..361 


485.000 


0 . 0 


0.0 


903.681 


0.362 


490.-000 


0.0 


0.0 


912.997 


2.536 


500.000 


0.0 


0.0 


931.630 


2 . 542 


510.000 


0 . 0 


0.0 


950.263 


2.548 


515.000 


0.0 


0.0 


959. 579 


12.388 


530.000 


0.0 


0.0 


987.528 


12. .427 


545.000 


0 . 0 


0.0 


1015.477 


12. 464 


555.000 


0.0 


Q.O 


1034. 109 


6 .- 244 -- 


575.000 


0.0 


0.0 


1071.374 


6.266 


595.000 


0 . 0 


0.0 


1108.640 


6.286 


600.000 


0.0 


0.0 


1117.956 


L..480 


635.000 


0.0 


0.0 


1183.170 


1.487 


660.000 


0.0 


0.0 


1229. 752 


1 .4-92 


670.000 


0.0 


G .0 


1248.384 


0 . 373 


705.000 


0.0 


0.0 


1313.598 


0.37 5 


730.000 


0.0 


0.0 


1360 . 180 


0 . 376 


740.000 


0.0 


0.0 


1378.812 


0.376 


775.000 


0 . 0 ■ 


0.0 


1444. 026 


0.377 


800.000 


0.0 


0.0 


1490.608 


0 . 378 . 


810.000 


0.0 


0.0 


1509.241 


0.0 



XLO 


RLT1 


PH 1 1 


RLT2 


PH 1.21 


400.000 


0.0 


0.0 


Q.O 


0.-0 


465.000 


0.0 


0.0 


Q.O 


0.0 


475.000 


0.0 


0.0 


Q.O 


oi. o: 


485.000 


0.0 


0.0 


0.0 


0.0 


490. 000 


0 . 0 


0.0 


0.0 


0.0 


500.000 


187.500 


1.082 


0.0 


Q.O 


510. COO 


382.778 


1.082 


0.0 


0.0 


515.000 


433.688 


5.256 


0.0 


Q.O 


530.000 


544.064 


5.256 


0.0 


0..0 


545.000 


625.860 


5.256 


0.0 


0.0 


555.000 


672. 289 


2.628 


0 . 0 


0 . 0 


575.000 


753.476 


2.628 


0.0 


0 . 0 


595. 000 


324. 555 


2.628 


0.0 


0.0 


600.000 


841.200 


0.618 


0.0 


0.0 


635.000 


948.366 


0.618 


0.0 


0.0 


660.000 


1013.887 


0.618 


0.0 


0.0 


670.000 


1045.730 


0 . 155 


0.0 


0 . 0 


705 .000 


1135.654 


0 . 155 


0.0 


0.0 


730.000 


1196. 856 


0.155 


Q.O 


0.0 


740.000 


1220.771 


0.155 


0.0 


0.0 


775 .000 


1302.382 


0.155 


0.0 


0.0 


800.000 


1359.016 


0.155 


Q.O 


0.0 


810. OCO 


1381.344 


0.0 


Q.O 


0.0 
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HEIGHT 


OFF THE 


BOTTOM = 


1197.000 


METERS 


XLO 


RADIUS 


PHT 


RLT3 


PH 1 3 


400. 300 


0. 0 


0.0 


1098.289 


0.0 


465.000 


0.0 


0.0 


1276. 761 


0. .196 


475.000 


0.0 


0 .0 


1304.219 


0..197 


485.000 


0.0 


Q.O 


1331.676 


Q. .19 7 


490.000 


0.0 


0.0 


1345.404 


1.3 34 


500.000 


0.0 


a -a 


1372. 862 


1.389 


510.000 


0. 0 


0.0 


1400. 319 


1.393 


515.000 


0.0 


0.0 


1414.048 


6. .775 


530.000 


0.0 


Q.O 


1455.233 


6.303: 


545.000 


0. 0 


q.q 


1496.419 


6.829 


555.000 


0.0 


0.0 


1523. 876 


3.423 


575.000 


0.0 


0.0 


1578.791 


2.43 8 


595.000 


0.0 


Q.O 


1633.705 


3.453 


600.000 


0.0 


Q.O 


1647.434 


Q..813 


635. 000 


0.0 


Q.O 


1743.534 


0..318 


660.000 


0. 0 


0.0 


1812. 177 


0.822 


670.000 


0.0 


0.0 


1839.635 


0. 206 


705. 000 


0.0 


0 .0 


1935.735 


0.207 


730.000 


0.0 


0.0 


2004.378 


0.207 


740.000 


0.0 


0.0 


2031. 835 


0.203 


775.000 


3. 0 


Q.O 


2127. 93 6 


0.203 


800.000 


0.0 


0.0 


2196. 579 


0.209 


810.000 


0.0 


Q.Q 


2224.036 


0.0. 



XLO 


RLT 1 


PHIL 


RL.T2 


PH 12 


400.000 


0.0 


0.0 


Q..0 


o:..o 


465.000 


0. 3 


Q.Q 


a.,o 


a.o 


475.000 


0.0 


Q.O 


Q.Q 


or. o: 


435 .000 


0.0 


Q.O 


0.0 


0.0 


490.000 


0 . 0 


0.0 


0.0 


0.0 


500.000 


0.0 


Q.O 


0..0 


0.0 


510.000 


527.833 


0.338 


0.0 


0.0 


515.000 


611.460 


1.643 


0.0 


0.0 


530.000 


784.567 


1.643 


0.0 


0.0 


545.000 


909.296 


1 .643 


0.0 


Q.O 


555.000 


979. 561 


0.821 


0 . 0 


0.0 


575.000 


1 101 .499 


0.821 


0 . 0 


0 . 0 


595.000 


1207.712 


0.821 


0.0 


0.0 


600.000 


1232.532 


0.193 


0.0 


0.0 


635.000 


1392.713 


0 . 193 


0.0 


0.0 


660. 000 


1496.637 


0.193 


0.0 


0.0 


670.000 


1536.438 


0.048 


0 . 0 


0.0 


705 .000 


1669.648 


0.048 


0.0 


0.0 


730.000 


1760. 222 


0.048 


0 . 0 


0.0 


740.000 


1795.598 


0.048 


0.0 


0.0 


775.000 


1916.269 


0.048 


0 . Q 


Q.Q 


800.000 


1999.965 


0.048 


0.0 


0.0 


810.000 


2032.955 


Q.O 


0.0 


0.0 
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COMPUTER PROGRAM 



IMPLICIT REAL*0 (A-H,0-Z) 

DIMENSION Y ( 2 » 5 ) . THTA ( 4) , A ( 4 ,4 1 , XLO (25) , ZPRT (20) , 
ISZ125, 10) ,SXL0{ 25) , SY2 ( 25 , 1 0 ) ,S PH I { 25 , 1 0 ) , 

2SRLTI (25,10) , S PH 1 1 ( 25 ,10 ), S RLT2 ( 25 , ID) , SPH 1 2 ( 25 , 1 0 ) , 
4SRLT3 (25,10) , SPH I 3 ( 25 , 1 J ) , FY (2 , 5 ) , AZDEL (1) 

FUZZ, Yi, Y2, DVT) = ( ( ( 2*GAM/ ( Y2* Y2 ) )*( Y2-< 3*RTDEL/ Y1 ) 
1*IGAM/Y2-ZZ) )/ (3~I D-ZZ ) +2* GAM/ Y2 ) ) + < 1 +-( 3 *RTDEL/ Y I ) 
2*tG AM/ { Y2*Y2)+< GAM/Y2-ZZJ -DVT/YL) ) ) /(3*(D-Z2)+ 
32*GAM/Y2) 

F2 122 , Yi , Y2 ) = ( Y2*Y2-(.3*RTD EL* ( GAM-( ZZ*Y2) )/YlJ )/( 3*Y2 
l*( D-ZZU2=GAM) 

THIS PROGRAM USES THE OPTIMAL RUNGE - KUTTA- INTEGRATION 
METHOD. 

READ ( 5 , 101) GRAV, BETA ,ALPH,.D r RTDEL, GAM, XAC 
101 FORMAT C8E1Q. 3 ) 

READ (5, 101) ZMA X , Z 0 

ZO IS THE DEPTH AT WhTCH THE SOURCE IS LOCATED. 

ZM A X IS THE SURFACE DEPTH— HERE IT IS ZERO. . 

READ (5,101) I THTA { J ) , J=1 ,4) 

READ ( 5, 101) ( ( A ( I , J ) , J = 1 ,4 ) ,1=1., 4) 

THE * A * AND «THTA* ARRAYS ARE USED FOR ; THEE RUNGEr 
KUTTA INTEGRATION. 

WRITE! 6 ,666 ) 

666 FORMAT Cl*) 

WRITE! 6,630) GRAV, BETA ,-ALPH , Q,,RTDE'L, GAM , XACE 

630 FORMAT ( .» 0 » , 7 ( 5X , El 0 .3 ). ), 

WR I T E 1 6 , 633 ) ZMAX,Z0 

633 FORMAT ( 1 5X , • ZMAX = ', E10.3 ,' ZO — '',E10.3.) 

WR IT E » 6 , 634 ) (THTA(J) »J = L,4) 

634 FO RH A T 1 4 ( 5 X , E 1 0 . 3 ) ) 

WRITE! 6, 635) ( ( A ( I ,J), J = I,-4) ,i:=L,4) 

635 FORMAT (4615.3 ) 

READ 15 , 130) NZDEL, -NSKP ,NXL0 ,NZPRT 

NZDEL IS THE NUMBER CF INCREMENTS FOR- I NGEGR ATT ON . . 

NS KP IS A DEGUGGING V ARI.ABL E USED TO LIMIT THE 
LINES OF OUTFUT. 

NX LO IS THE NUMBER OF BUBBLES OF' DIFFERENT INITIAL- 
RADII TO BE ANALYZED., 

NZPRT IS THE NUMBER CF DEPTHS FOR WH I CH- COMPAR I SON 
DATA IS TO BE OUTPUTTED. 

130 FORMAT C31 10) 

RE AD (5, 101) (A ZD EL ( J ) ,J = 1, NZDEL). 

AZDEL IS THE INCREMENT FOR INTEGRATION . , 

WRITE (6,800) NZDEL , NS KP, NXLO , NZPRT 
800 FORM A T ( L5X , * NZDEL = • ,1.8, • NSKP = • , 1 .8 , > * NXLO =- • , 
118,' NZPRT = ' , 18 ) 

WR.ITEI 6,802) (AZDEL(J) ,J = 1, NZDEL) 

802 FORMAT (' THE AZDEL ARE ' ,/,3El 5.3) 

READ (5,101) (XLO ( J ) , J = 1 , NXLO ) 

XLO IS THE INITIAL BUBBLE RADIUS. 

WRITE (6, 631) (XLG(J) ,-J = L,NXLO) 

631 FORMAT CO', 'THE INITIAL RADIUS AR E • , / , 5X , 10 El 2 . 3 ) 

READ 15, 101) ( Z P R T ( J ) , J = 1 ,N ZPRT ) 

ZPRT IS THE DEPTH FOR WHICH COMPARISON DATA IS TO BE 
OUTPUTTED. 

WRITE! 6,632) ( ZPRT ( J ) ,J=1, NZPRT) 

632 FORMAT CO', 'THE INPUT HEIGHTS ARE • , / , 5X, IDE 12. 3) 
GB=GRAV/8ETA 

C NOW COMPUTE THE INITIAL BUBBLE VELOCITY. 

DC 80 K=1 , NXLO 
D02=DSQRT ( XL 0 ( K) ) 

DD3 = RTDEL/DSGRT ( GE/XAC ) 

CALL PHIS (XLO ( K ) ,PHIC) 

VT 10=GB-XL0( K) “XLO( K) /ALPH 
VT2Q=DSQRT(G3*XL0(K)/XAC) 

VTO= ( (-1.0+(1.0+(4.0*VTLO*VT10 )/ ( VT20*VT20) ) 
1**0.500000) / 

2t C2.0*VT 10*VT10 ) / (,VT 2Q*VT20.) ) ) *VT1.0 
VO=VTO 

C GB IS GRAV DIVIDED BY BETA. 

SXLO ( K )=XLO( K ) * . 1D7 

C ALL VARIABLES BEGINNING WITH S ARE STORAGE ARRAYS 
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FOR THE FINAL OUTPUT. 

NOTE THAT RAO 1 1 ARE READ LN AS METERS BUT" OUTPUTTED 
IN MICRONS. 

y ( I, n = 0 . 

Y(2, I )=XL01K.) 

I F LG = 0 
DO 45 1=1,2 
DO 46 J = l,4 
FY!l,J)=0. 

46 CONTINUE 
45 CONTINUE 

WRITE! 6,636} GB , VTO r VG , VT1 0 , VT2.0 

636 FORMAT l * 1 * ,2X , ' G3 = *>EL2.5,' VTO. = ,ET2^ 5., f VO . == • 
1EI2.5, 

• 2' VT1Q = * ,E12* 5 , r VT20 = • , E1.2..5’), 

WRIT E ( 6, 63? ) (Y(J, L),J=I,2) 

637 FORMAT { 10X ,6 13X , E1Z« 5 ) ) 

Z = Z 0 

KK=i 
KKK=0 
HHT =0.0 

HHT IS THE HEIGHT OFF THE BOTTOM IT IS EQUAL. TO Z--Z0 

I S0=0 
VT= VTO 
12=0 

ZT LD Y= ZO 

WRITE! 6,600) VO, SXLOCK) , ZO , ZTLDY , VTO 
600 FORMAT! 'O' ,/////, 20X,' VO = ',EI0.3,' METERS/SEC ', 
120X , 

2' LO = * ,E10.3, r MICRONS ",/.,20X,'’ ZO =",F10.3, 

3* METERS *,/, 

420X, * ZTLDY = ' ,FIG.3,«' METERS ,! ,//,20X,' VO: TERMINAL 

5 ', 

6 E 1 0 . 3 , ' METERS /SEC* ) 

WRITE! 6,666) 

II = 0 
IS=0 

XXL0=XL0(K}*.LD7 
701 11=11+1 

ZDEL = AZDEL 1 1 } 

IS=I S+l 

IFdFLG.EQ.il GO TO 661 
99 DO 23 J= 1 , 4 

VT1=G3*Y!2,J)*Y(2, J) /ALPH 
VT2=DSQRTIGS*Y(2,J} /XAC) 

VT R2 =V T 1 * V T 1 / { VT 2- V T 2 ) 

VT=( (-1.+DSQRT! l.+4.*VTR2) ) / ( 2 . * VTR2 ) ) *VT1 
DVT IS DVT/DL !A DERIVATIVE). 

DVT = ( 3 *V Ti /Y (2 , J ) } / ( 1 + 4-VTR2 ) . 500000-VT / Y { 2 , J ) 

THE INTEGRATION IS DONE HERE. 

FY( 1 , J )=F1 ( Z+THTA! J ) *Z DE L , VT , Y ( 2 , J ) ,DVT) 

FY(2 , J )=F2 CZ+THTAt J } - ZDEL , VT , Y ( 2 , J ) } 

DO 2 5 N = 1 , 2 

Y(M, J + 1)=Y{N,1} + (A(J, 1)*FY(N,1)+A( J, 2) *FY( N ,2 ) +A! J,3) 
1*FY l N , 3 ) + 

2 A f J , 4 } -'- F Y ( N , 4 ) )*ZCEL 

WHEN THE RADIUS Y(2,J) IS LESS THAN ZERO, I FLG IS 
SET TO ONE AND THE INTEGRATING STEPS ARE SKIPPED 
FOR THE REST OF THE PROGRAM. 

I F { Y { N , J + 1 ) . LT .0.) GO TO 661 
25 CONTINUE 
23 CONTINUE 

Y 1 1 , 1 } =Y (1,5) 

Y( 2, 1 } =Y ( 2 , 5 ) 

GO TO 662 

661 Y ( 1 , 1 ) =1 . 

Yl 2, 11=0. 

I FLG= 1 

PH 1 = 0. 

662 Z=Z+ ZD EL 

IFtZ.GT . .1000-01 ) GO TO 999 
HHT =HHT+ ZDEL 
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I F ( IS.EQ.NSKP) GO TO 23 
GO TO 701 
28 IS=0 

I F l IFLG.EQ.l) GO TO 663 
VT1=G0*Y(2,1)*Y(2, 1)/ALPH 
VT2=DSQRT(G8*Y (2 ,1 J/XAC) 

VTR2 = VT 1*VT1 / ( VT2* VT2 ) 

VT = { (-1.+DSQRT ( 1 .+4.*VTR2) ) / ( 2 . *VTR2 ) ) *VT1 
GAM77=DEXP(Y(1 , 1 ) ) 

PH 1 = VT 0*PH 1 0/ ( VT*G AM77 ) 

C PHI IS THE GENERAL CASE DISTRIBUTION. 

663 RZ= ( D-Z) / ( D-ZO ) 

DO 1 1 = ( D-ZO)/(D-Z) 

DDl=DSQRT( { D— ZO ) / ( D- Z ) ) 

RLT1 = { DD11*XL0( K)*XLO< K)*XL0(K)+1. 5*ALPH*RTDEL 
1*IZ*Z-Z0*Z0)/(GB*(D-Z) ) ) 

IFIRLTl .LT..1D-18) RLT1=0. 

RLT1=RLT 1**0. 333 33 3333333 
PHI1=PHI0*RZ 

C PH 1 1 IS FOR THE VISCOUS DRAG MODEL. 

IF(RLTl.EQ.O.O) PHI1=0.0 

RLT2=( 001*002*002*002+003* ( (D-Z) *( 1.-001*001*001 )-3 . 
1*0*11.-001) ) ) 

IF(RLT2.LT. . 10-09) PLT2=0.. 

RLT2=RLT 2**0. 666666666666 
PH I 2 = PH I 0* DSQRT (RZ) 

PH 1 2 IS FOR THE TURBULENT DRAG, SURFACE ACTIVE 
MATERIALS MODEL. 

IF (RLT2. EQ. 0.0) PHI 2 = 0.0 
RZ3=RZ**0. 333333333333 
RLT3=XL0(K)/RZ3 
VT13=GB*RLT3*RLT3/ALPH 
VT23=DSQRT(GB*RLT3/XAC ) 

VTR3=VT13*VT13/ ( VT23*VT23) 

VT3=( (-1.+DSQRTI l.+A.*VTR3) )/(2.*VTR3) )*VT13 
PHI3=PHI0*VT0*RZ3/VT3 
C PH I 3 IS FOR THE NO GAS DIFFUSION CASE. 

Y2=Y( 2,1) *1.06 
XXXX=DABS ( HHT-ZPRT ( KK ) ) 

RRLT1=RLT1*. 107 
RRLT 2=RLT 2*. 107 
RRLT 3 = RLT 3*. 107 

IF(XXXX.GT.(0.1*ZDEL) ) GO TO 223 

KKK=KKK+ 1 

SZ(K,KKK)=Z 

SY2 ( K , KKK ) =Y ( 2 , 1 ) * . 1 07 
SPHI { K , KKK ) = PH I 
SRLT1 ( K,KKK)=RLT1* .107 
SPHI 1 ( K , KKK) =PH I 1 
$RLT2(K,KKK)=RLT2*.1D7 
SPH I 2 ( K , KKK ) =P H 1 2 
SRLT3 ( K , KKK) =R LT3* .107 
SPHI3(K,KKK)=PHI3 
KK=KK+ 1 
223 GO TO 701 
999 CONTINUE 
80 CONTINUE 

DO 81 KKK=1 , NZPRT 
WR I T E ( 6 , 82 ) 

82 FORMAT Cl') 

WR ITE ( 6, 83 ) ZPRT(KKK) 

83 FORMAT ( 'O’ ,/////, 30X, ' HEIGHT OFF THE BOTTOM = ', 
1F13. 3, ' METERS • ,///) ' 

WRITE ( 6,90) 

90 FORMAT ( 23X , ' XLO ',' RADIUS',' PHI ', 

1' P.LT3 ' , ' PH I 3 ' ) 

DO 85 K= 1 , NXL 0 

WRITE (6,86 ) SXLO ( K ) ,SY2 (K, KKK) , SPHI ( K , KKK ) , 

1SRLT3 ( K ,KKK) , SPH I 3 ( K , KKK ) 

85 CONTINUE 
WRITE (6 ,92 ) 

92 FORMAT ( • O' ,////) 
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WR I T E ( 6 , 9 1 ) 

91 FORMAT ( 23X » ' XLO RLT1 PHIL ' 

1* RLT2',« PHI 2') 

DO 87 K = 1 1 NXLO 

WRITE! 6 ,86) SXLC(K) ,SRLT1 ( K , KKK ) , SPH 1 1 { K , KKK ) , 
1SRLT2(K,KKK) ,SPHI2(K,KKK) 

86 FORMAT (23X,5FU. 3) 

87 CONTINUE 
81 CONTINUE 

STOP 

END 



SUBROUTINE PHIS(R,PHTO) 

IMPLICIT REAL*8 (A-H,0-Z) 

R IS MEASURED IN METERS 
IF(R .LT. .46500-03) GO TO 1 
I F ( (R.GE. .46 500-03) .AMD. (R.LT. .49000-03) ) 
IF( (R.GE. .49000-03) .AMD. (R.LT. . 5150D-03) ) 

I F ( (R.GE. .5150 D— 03 ) .AND. (R.LT. .5 55 00D-03) ) 
I F ( ( R.GE. .555 00-03 ) . AND. ( R . LT . .60000-03 ) ) 

I F ( (R.GE. .60000-03). AND. (R.LT. .67000-03) ) 

I F ( (R.GE. .67000-03). AND. ( R . LT .. 74000-03 ) ) 

I F ( ( R.GE. .74000-03 ) .AND. ( R .LT. . 8100D-03) ) 
IF ( R .G E . . 8 100 D— 03 ) GO TO 9 

1 PH I 0=0 . 

GO TO 10 

2 PH I 0 =1 . 

GO TO 10 

3 PH I 0=7 . 

GO TO 10 

4 PHI 0 = 34. 

GO TO 10 

5 PHI 0=17. 

GO TO 10 

6 PHI0=4 * 

GO TO 10 

7 PHIO= 1 . 

GO TO 10 

8 PH I 0=1. . 

GO TO 10 

9 PH 10 = 0 . 

10 CONTINUE 

RETURN 

END 



GO TO 2 
GO TO 3 
GO TO 4 
GO TO 5 
GO TO 6 
GO TO 7 
GO TO 8 
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